knitr::opts_chunk$set(
collapse = TRUE,
comment = "#>"
)library(Giotto)To run this vignette you need to install all the necessary Python modules.
This can be done manually, see https://rubd.github.io/Giotto_site/articles/installation_issues.html#python-manual-installation
This can be done within R using our installation tools (installGiottoEnvironment), see https://rubd.github.io/Giotto_site/articles/tut0_giotto_environment.html for more information.
# to automatically save figures in save_dir set save_plot to TRUE
temp_dir = '../results/' #'~/Temp/'
myinstructions = createGiottoInstructions(save_dir = temp_dir,
save_plot = FALSE,
show_plot = F)
#>
#> no external python path or giotto environment was specified, will check if a default python path is available
#>
#> A default python path was found: /opt/conda/bin/python3 and will be used
#>
#> If this is not the correct python path, either
#>
#> 1. use installGiottoEnvironment() to install a local miniconda python environment along with required modules
#>
#> 2. provide an existing python path to python_path to use your own python path which has all modules installedminimum requirements:
- matrix with expression information (or path to)
- x,y(,z) coordinates for cells or spots (or path to)
# giotto object
expr_path = system.file("extdata", "starmap_expr.txt.gz", package = 'Giotto')
loc_path = system.file("extdata", "starmap_cell_loc.txt", package = 'Giotto')
starmap_mini <- createGiottoObject(raw_exprs = expr_path,
spatial_locs = loc_path,
instructions = myinstructions)
#> Consider to install these (optional) packages to run all possible Giotto commands for spatial analyses: SPARK multinet RTriangle FactoMiner
#> Giotto does not automatically install all these packages as they are not absolutely required and this reduces the number of dependenciesHow to work with Giotto instructions that are part of your Giotto object:
- show the instructions associated with your Giotto object with showGiottoInstructions
- change one or more instructions with changeGiottoInstructions
- replace all instructions at once with replaceGiottoInstructions
- read or get a specific giotto instruction with readGiottoInstructions
Of note, the python path can only be set once in an R session. See the reticulate package for more information.
# show instructions associated with giotto object (starmap_mini)
showGiottoInstructions(starmap_mini)
#> $python_path
#> [1] "/opt/conda/bin/python3"
#>
#> $show_plot
#> [1] FALSE
#>
#> $return_plot
#> [1] TRUE
#>
#> $save_plot
#> [1] FALSE
#>
#> $save_dir
#> [1] "../results/"
#>
#> $plot_format
#> [1] "png"
#>
#> $dpi
#> [1] 300
#>
#> $units
#> [1] "in"
#>
#> $height
#> [1] 9
#>
#> $width
#> [1] 9
#>
#> $is_docker
#> [1] FALSEfilterDistributions(starmap_mini, detection = 'genes')filterDistributions(starmap_mini, detection = 'cells')filterCombinations(starmap_mini,
expression_thresholds = c(1),
gene_det_in_min_cells = c(50, 100, 200),
min_det_genes_per_cell = c(20, 28, 28))#> $results
#> threshold gene_detected_in_min_cells min_detected_genes_per_cell combination
#> 1: 1 50 20 50-20
#> 2: 1 100 28 100-28
#> 3: 1 200 28 200-28
#> removed_genes removed_cells
#> 1: 0 95
#> 2: 0 586
#> 3: 0 586
#>
#> $ggplot
starmap_mini <- filterGiotto(gobject = starmap_mini,
expression_threshold = 1,
gene_det_in_min_cells = 50,
min_det_genes_per_cell = 20,
expression_values = c('raw'),
verbose = T)
#> Number of cells removed: 95 out of 4000
#> Number of genes removed: 0 out of 28
starmap_mini <- normalizeGiotto(gobject = starmap_mini, scalefactor = 6000, verbose = T)
#>
#> first scale genes and then cells
starmap_mini <- addStatistics(gobject = starmap_mini)starmap_mini <- runPCA(gobject = starmap_mini, method = 'factominer')
#> hvg was not found in the gene metadata information, all genes will be used
#> Warning in runPCA_factominer(x = t_giotto(expr_values), scale = scale_unit, :
#> ncp > ncol(x), will be set to ncol(x)
screePlot(starmap_mini, ncp = 30)
#> PCA with name: pca already exists and will be used for the screeplotplotPCA(gobject = starmap_mini)
# 2D umap
starmap_mini <- runUMAP(starmap_mini, dimensions_to_use = 1:8)
plotUMAP(gobject = starmap_mini)
# 3D umap
starmap_mini <- runUMAP(starmap_mini, dimensions_to_use = 1:8, name = '3D_umap', n_components = 3)
plotUMAP_3D(gobject = starmap_mini, dim_reduction_name = '3D_umap')
#> create 3D plot
# 2D tsne
starmap_mini <- runtSNE(starmap_mini, dimensions_to_use = 1:8)
plotTSNE(gobject = starmap_mini)starmap_mini <- createNearestNetwork(gobject = starmap_mini, dimensions_to_use = 1:8, k = 25)
starmap_mini <- doLeidenCluster(gobject = starmap_mini, resolution = 0.5, n_iterations = 1000)
# 2D umap
plotUMAP(gobject = starmap_mini, cell_color = 'leiden_clus', show_NN_network = T, point_size = 2.5)
# 3D umap
plotUMAP_3D(gobject = starmap_mini, dim_reduction_name = '3D_umap', cell_color = 'leiden_clus')
#> create 3D plot
#> center label is not clear to see in 3D plot
#> You can shut it down with show_center_label = F
# 2D umap + coordinates
spatDimPlot(gobject = starmap_mini, cell_color = 'leiden_clus',
dim_point_size = 2, spat_point_size = 2.5)
# 3D umap + coordinates
spatDimPlot3D(gobject = starmap_mini, cell_color = 'leiden_clus', dim_reduction_name = '3D_umap')
#> center label is not clear to see in 3D plot
#> You can shut it down with show_center_label = F
#> Warning: 'layout' objects don't have these attributes: 'NA'
#> Valid attributes include:
#> 'font', 'title', 'uniformtext', 'autosize', 'width', 'height', 'margin', 'paper_bgcolor', 'plot_bgcolor', 'separators', 'hidesources', 'showlegend', 'colorway', 'datarevision', 'uirevision', 'editrevision', 'selectionrevision', 'template', 'modebar', 'meta', 'transition', '_deprecated', 'clickmode', 'dragmode', 'hovermode', 'hoverdistance', 'spikedistance', 'hoverlabel', 'selectdirection', 'grid', 'calendar', 'xaxis', 'yaxis', 'ternary', 'scene', 'geo', 'mapbox', 'polar', 'radialaxis', 'angularaxis', 'direction', 'orientation', 'editType', 'legend', 'annotations', 'shapes', 'images', 'updatemenus', 'sliders', 'colorscale', 'coloraxis', 'metasrc', 'barmode', 'bargap', 'mapType'
# heatmap and dendrogram
showClusterHeatmap(gobject = starmap_mini, cluster_column = 'leiden_clus')showClusterDendrogram(starmap_mini, h = 0.5, rotate = T, cluster_column = 'leiden_clus')gini_markers = findMarkers_one_vs_all(gobject = starmap_mini,
method = 'gini',
expression_values = 'normalized',
cluster_column = 'leiden_clus',
min_genes = 20,
min_expr_gini_score = 0.5,
min_det_gini_score = 0.5)
#>
#> start with cluster 1
#>
#> start with cluster 2
#>
#> start with cluster 3
#>
#> start with cluster 4
#>
#> start with cluster 5
#>
#> start with cluster 6
#>
#> start with cluster 7
#>
#> start with cluster 8
# get top 2 genes per cluster and visualize with violinplot
topgenes_gini = gini_markers[, head(.SD, 2), by = 'cluster']
violinPlot(starmap_mini, genes = topgenes_gini$genes, cluster_column = 'leiden_clus')
#> These genes have duplicates:
#> Egr1 Vip Vip Vip Pvalb Vip Sst Sst Egr2
# get top 6 genes per cluster and visualize with heatmap
topgenes_gini2 = gini_markers[, head(.SD, 6), by = 'cluster']
plotMetaDataHeatmap(starmap_mini, selected_genes = topgenes_gini2$genes,
metadata_cols = c('leiden_clus'))clusters_cell_types = c('cell A', 'cell B', 'cell C', 'cell D',
'cell E', 'cell F', 'cell G', 'cell H')
names(clusters_cell_types) = 1:8
starmap_mini = annotateGiotto(gobject = starmap_mini,
annotation_vector = clusters_cell_types,
cluster_column = 'leiden_clus',
name = 'cell_types')
# check new cell metadata
pDataDT(starmap_mini)
# visualize annotations
spatDimPlot(gobject = starmap_mini, cell_color = 'cell_types',
spat_point_size = 2, dim_point_size = 2)dimGenePlot3D(starmap_mini,
dim_reduction_name = '3D_umap',
expression_values = 'scaled',
genes = "Pcp4",
genes_high_color = 'red', genes_mid_color = 'white', genes_low_color = 'darkblue')
spatGenePlot3D(starmap_mini,
expression_values = 'scaled',
genes = "Pcp4",
show_other_cells = F,
genes_high_color = 'red', genes_mid_color = 'white', genes_low_color = 'darkblue')Create a grid based on defined stepsizes in the x,y(,z) axes.
starmap_mini <- createSpatialGrid(gobject = starmap_mini,
sdimx_stepsize = 200,
sdimy_stepsize = 200,
sdimz_stepsize = 20,
minimum_padding = 10)
showGrids(starmap_mini)
#> The following grids are available: spatial_grid
#> [1] "spatial_grid"
# visualize grid
spatPlot2D(gobject = starmap_mini, show_grid = T, point_size = 1.5)Only the method = delaunayn_geometry can make 3D Delaunay networks. This requires the package geometry to be installed.
plotStatDelaunayNetwork(gobject = starmap_mini, maximum_distance = 200,
method = 'delaunayn_geometry')
#> Warning: Ignoring unknown parameters: binwidth, bins, padstarmap_mini = createSpatialNetwork(gobject = starmap_mini, minimum_k = 2,
maximum_distance_delaunay = 200,
method = 'Delaunay',
delaunay_method = 'delaunayn_geometry')
starmap_mini = createSpatialNetwork(gobject = starmap_mini, minimum_k = 2,
method = 'kNN', k = 10)
showNetworks(starmap_mini)
#> The following images are available: Delaunay_network kNN_network
#> [1] "Delaunay_network" "kNN_network"
# visualize the two different spatial networks
spatPlot(gobject = starmap_mini, show_network = T,
network_color = 'blue', spatial_network_name = 'Delaunay_network',
point_size = 2.5, cell_color = 'leiden_clus')
spatPlot(gobject = starmap_mini, show_network = T,
network_color = 'blue', spatial_network_name = 'kNN_network',
point_size = 2.5, cell_color = 'leiden_clus')Identify spatial genes with 3 different methods:
- binSpect with kmeans binarization (default)
- binSpect with rank binarization
- silhouetteRank
Visualize top 4 genes per method.
km_spatialgenes = binSpect(starmap_mini)
#>
#> This is the single parameter version of binSpect
#> 1. matrix binarization complete
#>
#> 2. spatial enrichment test completed
#>
#> 3. (optional) average expression of high expressing cells calculated
#>
#> 4. (optional) number of high expressing cells calculated
spatGenePlot(starmap_mini, expression_values = 'scaled',
genes = km_spatialgenes[1:4]$genes,
point_shape = 'border', point_border_stroke = 0.1,
show_network = F, network_color = 'lightgrey', point_size = 2.5,
cow_n_col = 2)
rank_spatialgenes = binSpect(starmap_mini, bin_method = 'rank')
#>
#> This is the single parameter version of binSpect
#> 1. matrix binarization complete
#>
#> 2. spatial enrichment test completed
#>
#> 3. (optional) average expression of high expressing cells calculated
#>
#> 4. (optional) number of high expressing cells calculated
spatGenePlot(starmap_mini, expression_values = 'scaled',
genes = rank_spatialgenes[1:4]$genes,
point_shape = 'border', point_border_stroke = 0.1,
show_network = F, network_color = 'lightgrey', point_size = 2.5,
cow_n_col = 2)
silh_spatialgenes = silhouetteRank(gobject = starmap_mini) # TODO: suppress print output
spatGenePlot(starmap_mini, expression_values = 'scaled',
genes = silh_spatialgenes[1:4]$genes,
point_shape = 'border', point_border_stroke = 0.1,
show_network = F, network_color = 'lightgrey', point_size = 2.5,
cow_n_col = 2)Identify robust spatial co-expression patterns using the spatial network or grid and a subset of individual spatial genes.
1. calculate spatial correlation scores
2. cluster correlation scores
# 1. calculate spatial correlation scores
ext_spatial_genes = km_spatialgenes[1:20]$genes
spat_cor_netw_DT = detectSpatialCorGenes(starmap_mini,
method = 'network',
spatial_network_name = 'Delaunay_network',
subset_genes = ext_spatial_genes)
# 2. cluster correlation scores
spat_cor_netw_DT = clusterSpatialCorGenes(spat_cor_netw_DT,
name = 'spat_netw_clus', k = 6)
heatmSpatialCorGenes(starmap_mini, spatCorObject = spat_cor_netw_DT,
use_clus_name = 'spat_netw_clus')
netw_ranks = rankSpatialCorGroups(starmap_mini,
spatCorObject = spat_cor_netw_DT,
use_clus_name = 'spat_netw_clus')
top_netw_spat_cluster = showSpatialCorGenes(spat_cor_netw_DT,
use_clus_name = 'spat_netw_clus',
selected_clusters = 6,
show_top_genes = 1)
cluster_genes_DT = showSpatialCorGenes(spat_cor_netw_DT,
use_clus_name = 'spat_netw_clus',
show_top_genes = 1)
cluster_genes = cluster_genes_DT$clus; names(cluster_genes) = cluster_genes_DT$gene_ID
starmap_mini = createMetagenes(starmap_mini, gene_clusters = cluster_genes, name = 'cluster_metagene')
spatCellPlot(starmap_mini,
spat_enr_names = 'cluster_metagene',
cell_annotation_values = netw_ranks$clusters,
point_size = 1.5, cow_n_col = 3)hmrf_folder = paste0(temp_dir,'/','11_HMRF/')
if(!file.exists(hmrf_folder)) dir.create(hmrf_folder, recursive = T)
# perform hmrf
my_spatial_genes = km_spatialgenes[1:20]$genes
HMRF_spatial_genes = doHMRF(gobject = starmap_mini,
expression_values = 'scaled',
spatial_genes = my_spatial_genes,
spatial_network_name = 'Delaunay_network',
k = 6,
betas = c(10,2,2),
output_folder = paste0(hmrf_folder, '/', 'Spatial_genes/SG_top20_k6_scaled'))
#>
#> expression_matrix.txt already exists at this location, will be overwritten
#>
#> spatial_genes.txt already exists at this location, will be overwritten
#>
#> spatial_network.txt already exists at this location, will be overwritten
#>
#> spatial_cell_locations.txt already exists at this location, will be overwritten
#> [1] "/opt/conda/bin/python3 /usr/local/lib/R/site-library/Giotto/python/reader2.py -l \"../results//11_HMRF//Spatial_genes/SG_top20_k6_scaled/spatial_cell_locations.txt\" -g \"../results//11_HMRF//Spatial_genes/SG_top20_k6_scaled/spatial_genes.txt\" -n \"../results//11_HMRF//Spatial_genes/SG_top20_k6_scaled/spatial_network.txt\" -e \"../results//11_HMRF//Spatial_genes/SG_top20_k6_scaled/expression_matrix.txt\" -o \"../results//11_HMRF//Spatial_genes/SG_top20_k6_scaled/result.spatial.zscore\" -a test -k 6 -b 10 2 2 -t 1e-10 -z none -s 100 -i 100"
# check and select hmrf
for(i in seq(10, 14, by = 2)) {
viewHMRFresults2D(gobject = starmap_mini,
HMRFoutput = HMRF_spatial_genes,
k = 6, betas_to_view = i,
point_size = 2)
}
#> [1] "/opt/conda/bin/python3 /usr/local/lib/R/site-library/Giotto/python/get_result2.py -r \"../results//11_HMRF//Spatial_genes/SG_top20_k6_scaled/result.spatial.zscore\" -a test -k 6 -b 10"
#> Warning in system(command = result_command, intern = T): running command '/opt/
#> conda/bin/python3 /usr/local/lib/R/site-library/Giotto/python/get_result2.py -r
#> "../results//11_HMRF//Spatial_genes/SG_top20_k6_scaled/result.spatial.zscore" -a
#> test -k 6 -b 10' had status 1
#> [1] "/opt/conda/bin/python3 /usr/local/lib/R/site-library/Giotto/python/get_result2.py -r \"../results//11_HMRF//Spatial_genes/SG_top20_k6_scaled/result.spatial.zscore\" -a test -k 6 -b 12"
#> Warning in system(command = result_command, intern = T): running command '/opt/
#> conda/bin/python3 /usr/local/lib/R/site-library/Giotto/python/get_result2.py -r
#> "../results//11_HMRF//Spatial_genes/SG_top20_k6_scaled/result.spatial.zscore" -a
#> test -k 6 -b 12' had status 1
starmap_mini = addHMRF(gobject = starmap_mini,
HMRFoutput = HMRF_spatial_genes,
k = 6, betas_to_add = c(12),
hmrf_name = 'HMRF')
#> [1] "/opt/conda/bin/python3 /usr/local/lib/R/site-library/Giotto/python/get_result2.py -r \"../results//11_HMRF//Spatial_genes/SG_top20_k6_scaled/result.spatial.zscore\" -a test -k 6 -b 12"
#> Warning in system(command = result_command, intern = T): running command '/opt/
#> conda/bin/python3 /usr/local/lib/R/site-library/Giotto/python/get_result2.py -r
#> "../results//11_HMRF//Spatial_genes/SG_top20_k6_scaled/result.spatial.zscore" -a
#> test -k 6 -b 12' had status 1
# visualize selected hmrf result
giotto_colors = Giotto:::getDistinctColors(6)
names(giotto_colors) = 1:6
spatPlot(gobject = starmap_mini, cell_color = 'HMRF_k6_b.12',
point_size = 3, coord_fix_ratio = 1, cell_color_code = giotto_colors)set.seed(seed = 2841)
cell_proximities = cellProximityEnrichment(gobject = starmap_mini,
cluster_column = 'cell_types',
spatial_network_name = 'Delaunay_network',
adjust_method = 'fdr',
number_of_simulations = 1000)
# barplot
cellProximityBarplot(gobject = starmap_mini,
CPscore = cell_proximities,
min_orig_ints = 5, min_sim_ints = 5)## heatmap
cellProximityHeatmap(gobject = starmap_mini, CPscore = cell_proximities,
order_cell_types = T, scale = T,
color_breaks = c(-1.5, 0, 1.5),
color_names = c('blue', 'white', 'red'))
# network
cellProximityNetwork(gobject = starmap_mini, CPscore = cell_proximities,
remove_self_edges = T, only_show_enrichment_edges = T)
# network with self-edges
cellProximityNetwork(gobject = starmap_mini, CPscore = cell_proximities,
remove_self_edges = F, self_loop_strength = 0.3,
only_show_enrichment_edges = F,
rescale_edge_weights = T,
node_size = 8,
edge_weight_range_depletion = c(1, 2),
edge_weight_range_enrichment = c(2,5))pDataDT(starmap_mini)
#> cell_ID nr_genes perc_genes total_expr leiden_clus cell_types
#> 1: cell_2944 28 100.00000 186.0057 6 cell F
#> 2: cell_17523 28 100.00000 179.2093 4 cell D
#> 3: cell_30080 27 96.42857 174.9019 1 cell A
#> 4: cell_941 28 100.00000 188.7800 6 cell F
#> 5: cell_28514 28 100.00000 173.0342 2 cell B
#> ---
#> 3901: cell_26011 28 100.00000 195.5315 4 cell D
#> 3902: cell_31539 28 100.00000 160.8629 1 cell A
#> 3903: cell_4185 28 100.00000 184.1559 2 cell B
#> 3904: cell_3851 28 100.00000 198.2840 5 cell E
#> 3905: cell_21337 28 100.00000 123.7624 3 cell C
#> HMRF_k6_b.12
#> 1: <NA>
#> 2: <NA>
#> 3: <NA>
#> 4: <NA>
#> 5: <NA>
#> ---
#> 3901: <NA>
#> 3902: <NA>
#> 3903: <NA>
#> 3904: <NA>
#> 3905: <NA>
# Option 1
spec_interaction = "cell D--cell H" # needs to be in alphabetic order! first D, then H
cellProximitySpatPlot2D(gobject = starmap_mini,
interaction_name = spec_interaction,
show_network = T,
cluster_column = 'cell_types',
cell_color = 'cell_types',
cell_color_code = c('cell H' = 'lightblue', 'cell D' = 'red'),
point_size_select = 4, point_size_other = 2)
# Option 2: create additional metadata
starmap_mini = addCellIntMetadata(starmap_mini,
spatial_network = 'Delaunay_network',
cluster_column = 'cell_types',
cell_interaction = spec_interaction,
name = 'D_H_interactions')
spatPlot(starmap_mini, cell_color = 'D_H_interactions', legend_symbol_size = 3,
select_cell_groups = c('other_cell D', 'other_cell H', 'select_cell D', 'select_cell H'))
# create cross section
starmap_mini = createCrossSection(starmap_mini,
method="equation",
equation=c(0,1,0,600),
extend_ratio = 0.6)
# show cross section
insertCrossSectionSpatPlot3D(starmap_mini, cell_color = 'leiden_clus',
axis_scale = 'cube',
point_size = 2)
#> create 3D plot
insertCrossSectionGenePlot3D(starmap_mini, expression_values = 'scaled',
axis_scale = "cube",
genes = "Slc17a7")
# for cell annotation
crossSectionPlot(starmap_mini,
point_size = 2, point_shape = "border",
cell_color = "leiden_clus")
crossSectionPlot3D(starmap_mini,
point_size = 2, cell_color = "leiden_clus",
axis_scale = "cube")
#> create 3D plot
# for gene expression
crossSectionGenePlot(starmap_mini,
genes = "Slc17a7",
point_size = 2,
point_shape = "border",
cow_n_col = 1.5,
expression_values = 'scaled')
crossSectionGenePlot3D(starmap_mini,
point_size = 2,
genes = c("Slc17a7"),
expression_values = 'scaled')viewer_folder = paste0(temp_dir, '/', 'Mouse_cortex_viewer')
# select annotations, reductions and expression values to view in Giotto Viewer
exportGiottoViewer(gobject = starmap_mini, output_directory = viewer_folder,
factor_annotations = c('cell_types',
'leiden_clus',
'HMRF_k6_b.12'),
numeric_annotations = 'total_expr',
dim_reductions = c('umap'),
dim_reduction_names = c('umap'),
expression_values = 'scaled',
expression_rounding = 3,
overwrite_dir = T)
#> output directory already exits, files will be overwritten
#> write cell and gene IDs
#> write physical centroid locations
#> write annotation data for: cell_types
#> write annotation data for: leiden_clus
#> write annotation data for: HMRF_k6_b.12
#> write annotation data for: total_expr
#> write annotation data for: umap for umap
#> write expression values
#> finished writing giotto viewer files to ../results//Mouse_cortex_viewer
#>
#> ================================================================
#> Next steps. Please manually run the following in a SHELL terminal:
#> ================================================================
#> cd ../results//Mouse_cortex_viewer
#> giotto_setup_image --require-stitch=n --image=n --image-multi-channel=n --segmentation=n --multi-fov=n --output-json=step1.json
#> smfish_step1_setup -c step1.json
#> giotto_setup_viewer --num-panel=2 --input-preprocess-json=step1.json --panel-1=PanelPhysicalSimple --panel-2=PanelTsne --output-json=step2.json --input-annotation-list=annotation_list.txt
#> smfish_read_config -c step2.json -o test.dec6.js -p test.dec6.html -q test.dec6.css
#> giotto_copy_js_css --output .
#> python3 -m http.server
#> ================================================================
#>
#> Finally, open your browser, navigate to http://localhost:8000/. Then click on the file test.dec6.html to see the viewer.
#>
#>
#> For more information, http://spatialgiotto.rc.fas.harvard.edu/giotto.viewer.setup3.html